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The rheology of 2D steady surface flow of cohesionless cylinders in a rotating drum is investigated 
through Non Smooth Contact Dynamics simulations. Profile of volume fraction, translational and 
angular velocity, rms velocity, strain rate and stress tensor were measured at the midpoint along the 
length of the surface flowing layer where the flow is generally considered as steady and homogeneous. 
Analysis of these data and their inter-relations suggest the local inertial number - defined as the 
ratio between local inertial forces and local confinement forces - to be the relevant dimensionless 
parameter to describe the transition from the quasi-static part of the packing to the flowing part 
at the surface of the heap. Variations of the components of the stress tensor as well as the ones of 
rms velocity as a function of the inertial number are analysed within both the quasi-static and the 
flowing phases. Their implications are discussed. 

PACS numbers: 45.70.-n, 83.70.Fn, 46.10,+z 



I. INTRODUCTION 

Granular media present numbers of interesting and un- 
usual behaviours: They can flow as liquids, but, under 
some circumstances, they can jam and resist to external 
shear stress without deforming. Understanding rheology 
of granular systems has thus developed along two major 
themes: The Rapid flows - gaseous-like - regime where 
grains interact through binary collisions, are generally 
described in the framework of the kinetic theory 0,0,13; 
The slow flow - solid-like - regime where grain inertia 
is negligible is most commonly described using the tools 
of soil mechanics and plasticity theory jj]. In between 
these two regimes there exists a dense flow - liquid- 
like - regime where grain inertia becomes important but 
contacts between grains are kept. Rheology of this last 
regime has been widely investigated experimentally, nu- 
merically and theoretically (see [f| for a review), but still 
remains far from being understood. Several models have 
been proposed recently to describe dense granular flows 
by accounting for non-local effects 0, 0,Hi M, E3- by 
adapting kinetic theory pHII^lp i , b y modelling dense 
flows as partially fluidized flows |l4lll5l] or by considering 
them as quasi-static flows where the mean motion results 
from transient fractures modelled as self activated pro- 
cess but, to our knowledge, none of them 
succeed to account for all the features experimentally ob- 
served. 

The most spectacular manifestation of this solid/liquid 
duality occurs during an avalanche when a thin layer of 
grains starts to roll at the surface of the packing, most of 
the grains remaining apparently static. The global evo- 
lution of such surface flows can be captured by mod- 
els derived from non-linear physics [2(], |^, ^ or fluid 
mechanics [53, 13 HE H3- However, some experimen- 



tal results remain unexplained: For instance, experimen- 
tal veloci ty prof iles measured in two-dimensional (2D) 
flo ws pl l28 l l2 H or three dimensional (3D) flows plIM 
I3H l32| clearly exhibit the selection of a constant velocity 
gradient within the flowing layer while momentum bal- 
ance implies that the shear stress increases linearly with 
depth. This observation is incompatible with any local 
and one-to-one stress/strain constitutive relations. Re- 
cent experiments |33l ] have provided evidence of 'jammed' 
aggregates embedded in the avalanche. These 'solid' clus- 
ters are found to be power-law distributed without any 
characteristic length-scales, and may explain the fail- 
ure of present models. But a clear understanding of the 
avalanche rheology is still missing. 

The purpose of this paper is to investigate the sur- 
face flows rheology through numerical simulations of 2D 
'minimal' granular systems made of cohesionless weakly 
polydisperse cylinders confined in a slowly rotating drum. 
Those allow us to track the evolution of quantities like 
stress that are not accessible in real experiments. More- 
over, they allow to get rid of artefacts such as the friction 
of beads on the lateral boundaries of the drum that may 
confuse the interpretation of an experiment. The numer- 
ical simulation were performed using contact dynamic 
methods 0, |3f| based on a fully implicit resolution of 
the contact forces, without any resort to regularization 
schemes. At a given step of evolution, all the kinematic 
constraints within the packing are simultaneously taken 
into account together with the equations of motions to 
determine all the contact forces in the packing. This al- 
lows to deal properly with nonlocal momentum transfers 
implied in multiple collisions, contrary to Molecular Dy- 
namics schemes traditionally used that reduce the system 
evolution to a succession of binary collisions. 

The simulation scheme and the description of the sim- 
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ulated systems are detailed in Sec. II. In Sec. Ill, we 
report comprehensive analysis of volume fraction and ve- 
locity (translational and rotational) profiles at the center 
of the drum. They are compared with experimental data 
available in the literature. Stress analysis and implica- 
tions on the rheology of free surface flows are discussed 
in Sec. IV. and Section V respectively. Sec. VI is focussed 
on the analysis of both the translational and rotational 
velocity fluctuations. Finally, Sec. VII summarizes our 
findings. 



II. SIMULATION METHODOLOGY 

For this study, we simulate granular systems similar 
to those investigated experimentally in [27l l28l I29I l33l| . 
We model a 2D rotating drum of diameter D equal 
to 450 mm half-filled with 7183 rigid disks of density 
po = 2.7 g.cm~ 3 and diameter uniformly distributed be- 
tween 3 mm and 3.6 mm. This weak polydispersity pre- 
vents any 2D ordering effect that may induce nongeneric 
effects. Normal restitution coefficient between two disks 
(respectively between disks and drum) is set to 0.46 (re- 
spectively 0.46) and the friction coefficient to 0.4 (respec- 
tively 0.95). Normal restitution coefficients and disk/disk 
friction coefficient were chosen to mimic the experimen- 
tal flows of aluminium beads investigated in p7l.l28| . The 
drum/disk friction coefficient was set close to 1 to prevent 
sliding at the drum boundary. 

Numerical simulations dedicated to evolution ofgran- 
ular media can be based either on explicit |36Ll37l l38. 39] 
or implicit 0H! 113 

method. One of the drawbacks of 
explicit models is to reduce non-local momentum trans- 
fers implied in multiple collisions to a succession of binary 
collisions. Moreover, numerical instabilities can occur in 
granular flows. They are corrected either by introducing 
some artificial viscosity or by reducing the size of the 
time step. The Non Smooth Contact Dynamics method 
used here is implicit. It provides a nonsmooth formula- 
tion of the bodys impenetrability condition, the collision 
rules and the dry Coulomb friction law. The method is 
extensively described in l4l| . and briefly explained below. 

Firstly, equations of motion are written for a collection 
of rigid bodies and discretized by a time integrator |42j |. 
The interaction problem is then solved at contact level 
(local level) rather than at particle level (global level) as 
commonly performed in explicit methods. In other words, 
equations are written in term of relative velocities u a and 
local impulsions r Q defined at each contact point indexed 
by a. The impenetrability condition evoked previously 
means that particles candidates for contact should not 
cross the boundaries of antagonist's bodies. We consider 
also that contacting bodies do not attract each other, i.e. 
that the reaction force is positive, and vanishes when the 
contact vanishes. This can be summarized in the follow- 
ing so-called velocity Signorini Condition: 

u„ > r n > u n .r n = 0, (1) 



where the index n denotes the normal component of the 
various quantities (index a is omitted). Let us note that 
this philosophy is different from what is used in explicit 
methods, where normal forces are usually proportional to 
the penetration between two particles. The dry frictional 
law is the Coulomb's one for which the basic features are: 
The friction force lies in the Coulomb's cone (| \r t \ | < fj,r n , 
(i friction coefficient), and if the sliding relative velocity is 
not equal to zero, its direction is opposed to the friction 
force {\\r t \\ — ^r n ). This summarized in the following 
relation: 

||n||<A""n IKH 7^ r t = -A""nTT-TT ( 2 ) 

Ikll 

For rigid bodies we also need to adopt a collision law 
because the velocity Signorini condition does not give 
enough information. We adopt the Newton restitution 
law, u+ = —e n u~, realistic for collection of disks. The 
reader can refer to |43| for more explanations about col- 
lision laws. Time discretization of equations of motion 
- where the global contact forces are the only missing 
quantities to determine the motion of each bead - lead to 
the following scheme: 

(Wr-u=b 

1 law a [u a , r a ] = .true., a = l,n c ^ ' 

where u and r denotes the vectors containing the rela- 
tive velocity and the mean contact impulse for all the 
contact points respectively. The matrix W is the Delas- 
sus operator 0] that contains all the local informations 
(local frames and contact points) as well as the infor- 
mations related to the contacts connectivity. The right 
hand side of first line in Eq. [31 represents the free rela- 
tive velocity calculated by only taking into account the 
external forces. The operator law a encodes the friction- 
contact law which should be satisfied by each component 
of couple (u Q ,r Q ); n c denotes the number of contact. 
Systems of Eq. can be solved by a classical non lin- 
ear Gau/3-Seidel algorithm or a Conjugate Projected 
Gradient one . This two algorithms benefit from par- 
allel versions poi |47| which show their efficiency in the 
simulations of large systems. Information from this local 
level, the contact level, is transfered to the global level, 
the grain level and the configuration of the system is up- 
dated. 

The procedure to achieve a numerical experiment is 
the following: All the disks are placed in an immobile 
drum; Once the packing is stabilized, a constant rota- 
tion speed £1 (ranging from 2 rpm to 15 rpm) is imposed 
to the drum; After one round, a steady continuous sur- 
face flow is reached (This has been checked by looking 
at the time evolution of the total kinetic energy within 
the packing over the next round); One starts then to 
capture 400 snapshots equally distributed over a rota- 
tion of the drum. The time-step is set to 6.10 -3 s. The 
number of time-steps necessary to achieve an experiment 
ranges from 4.10 3 to 10 4 depending on the rotating speed. 
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FIG. 1: Typical snapshot of the steady surface flows in the 
simulated 2D rotating drum. Its diameter and its rotation 
speed are respectively Do = 450 mm and £7 = 6 rpm. It is 
filled with 7183 disks of density po = 2.7 g.cm -3 and diame- 
ter uniformly distributed in the interval [3 mm; 3.6 mm], the 
disk/disk coefficient of restitution and disk/disk friction co- 
efficient are set respectively to e n = 0.46 and = 0.2. the 
disk/drum coefficient of restitution and disk/drum friction 
coefficient are set respectively to e„ = 0.46 and fi° = 0.95. 

All simulations have been performed with LMGC90 soft- 
ware . On a SGI Origin 3800 with 16 processors, about 
20 h are required to achieve one of these simulations. 

A typical snapshot of the simulated granular packing is 
shown in Fig. 2] For each bead of each of the 400 frames 
within a given numerical experiment, one records the po- 
sition x of its center of mass, the "instantaneous" veloc- 
ity c of this center of mass measured over a time win- 
dow St = 6.10 -3 s and its angular velocity w. Voronoi 
tessellation was used to define the local volume fraction 
associated to each bead (see e.g. 33]). The components 
of contact stress tensor a associated to each bead i arc 
computed as 0, |5(j : 

^^a^E^/*' a-- 9 6 {1,2}, (4) 

where Vi is the volume of the Voronoi polyedra associ- 
ated to the bead i, Fji the contact force between i and 
j, and Xjj = Xj — Xj. In all the following, these quantities 
are nondimcnsionalizcd: Calling g the gravity constant 
and d the mean disk diameter, distan ces, time, velocities 
and stresses are given in units of d, \J d/g, \J~gd and pogd 
respectively. In this paper, we concentrate on the contin- 
uum scale by looking at profiles of the time and space 
averaged quantities. Statistical analysis of these quanti- 
ties at the grain scale will be presented in a separated 
paper. 

In rotating drum geometries, the surface flow is not 



fully developed. The frame of study should now be cho- 
sen appropriately. One thus define the frame 3? rotating 
with the drum that coincides with the reference frame 
= (e K ,e 2 ) fixed in the laboratory, so that e x (resp. 
e z ) is parallel (resp. perpendicular) to the free surface 
(see Fig. j^Jj . In the frame 3?, the flow can be con- 
sidered as quasi-homogeneous at the center of the drum, 
e.g. within the elementary slice £ (see Fig. ^) 20 beads 
diameter wide, parallel to e z located at x = 0. This slice 
is divided into layers of one mean bead diameter wide 
parallel to the flow. The given value of a given contin- 
uum quantity a(z) (volume fraction, velocity, stress...) at 
depth z is then defined as the average of the correspond- 
ing quantity defined at the grain scale over all the beads 
in all the 400 frames of the sequence whose center of mass 
is inside the layer. 

III. KINEMATIC ANALYSIS 

A. Volume fraction profile 

Let us first focus on volume fraction profiles within 
the packing. Figure displays the volume fraction pro- 
file measured for = 6 rpm. To check the homogene- 
ity of the flow with regard to v, the elementary slice 
£ was translated of an increment of 5 bead diameters 
in both positive x and negative x. The volume fraction 
profile is found to be invariant under infinitesimal trans- 
lation along e x . At the free surface, v drops quickly to 
zero within a small zone of thickness around three/four 
beads diameter independent of £1. In all the following, 
the free surface boundary is set at the lower boundary 
of this small region (mixed line in Fig. [2J, defined at 
the point where v becomes larger than 0.7. At the drum 
boundary, v jumps also to a much smaller value within 
a small region about two/three beads diameters thick, 
which should be attributed to the presence of the smooth 
drum boundary. Apart fr21.8om these two narrow re- 
gions, the volume fraction v is almost constant within 
the drum, around the random close packing (RCP) value 
v rcp ^ Q_g2. A closer look (Inset of Fig. [2J suggests 
that v is constant within the static phase, and decreases 
weakly within the flowing layer as defined from the ve- 
locity profile in next section. Such behaviour is expected 
since granular systems should dilate before being able 
to deform. However, this decreasing is very small and 
compressible effects can thus be neglected with regard to 
momentum balance, even if they may significantly alter 
the local flow rheology [33| . 

B. Velocity profiles 

As expected for a quasi-homogeneous flow, the normal 
component v z of the velocity was found to be negligible 
compared to the tangential component v x at any depth 
z. Figure |31 depicts both the streamwise velocity profile 
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FIG. 2: Volume fraction profile viz) (averaged over 400 
frames) at the center of the drum obtained for Q — 6 rpm. 
The errorbars correspond to a 99% confident interval. The 
vertical mixed line locates the free surface boundary, - set 
to the point where v crosses the value 0.7 -, from which v 
drops quickly down to zero. The vertical dotted line locates 
the static/flowing interface, defined from the streamwise ve- 
locity profile (see Fig. |3J • Apart from these two zones, v ap- 
pears almost constant, around 0.8. Inset, zoom in the "con- 
stant region" enhancing the small variations of v within the 
two phases. 



0.25 
0.2 



N 0.15 

X 

r% 0.1 



0.05 




10° 










N 

ro 










- 10~ 1 


X 

> 


















it 




10- 3 












10~ 4 






Z I If 






-60 


-40 


-20 IT 














(b) 















-60 



-40 



-20 



v x (z) (Fig. |3t) an d the shear rate profile d z v x (z) (Fig. 

for £1 = 6 rpm. Both these profiles were found to be 
invariant under infinitesimal translation along e x . Two 
phases can clearly be observed: A flowing layer exhibiting 
a linear velocity profile and a static phase experiencing 
creep motion where both v x , and dv x /dz decay exponen- 
tially with depth (see Inset of Fig. 0b)- Such shapes are 
very similar to the ones observed experimentally in 2D 
flows IH E3 as well as in 3D flows jH |H |M H2 . 
The interface between the two phases can then be defined 
by extrapolating the linear velocity profile of the flowing 
phase to zero (see Fig. EK) - The flowing layer thickness 
H can then deduced. 

Velocity profiles measured for = 6 rpm at five dif- 
ferent locations x are represented in Fig. 0] At these lo- 
cations, d x H is no more equal to zero. The width of the 
elementary slice £ has thus been decreased to two beads 
diameters in order to minimize this drift. The shape of 
the velocity profile remains the same in these locations, 
with a clear linear velocity profile within the flowing layer 
and an exponentially decaying velocity within the static 
phase. Both the characteristic decay length A of the ex- 
ponential creep within the static phase and the constant 
velocity gradient 70 within the flowing layer are observed 
to be independent of the precise location x for a given 
value of fi. 

The "natural" control parameter in our experiment is 
the rotating speed Q. However, comparisons between ex- 
periments in heap geometry and rotating drum geometry 
.5] suggest that the main control parameter for the sur- 
face flow is the non-dimensionalized flow rate Q, defined 



FIG. 3: Profile of (a) streamwise velocity v x (z) and (b) 
streamwise velocity gradient d z v x (z) (averaged over 400 
frames) at the center of the drum as obtained for = 6 rpm. 
The errorbars correspond to a 99% confident interval. The 
profile of the velocity gradient (resp. of the velocity gradi- 
ent) is linear (resp. constant) within the flowing layer. The 
plain straight line in sub-figure (a) (resp. in sub-figure (b)) 
corresponds to a linear fit: v x = 7(2 + H) (resp. a constant 
d z v x = 7) where 7 ~ 0.15. The vertical mixed line locates the 
free surface boundary (see Fig.|2J. The flowing/static interface 
(dotted line) is defined from the depth where the straight line 
intersects the z-axis in sub- figure (a). The flowing layer thick- 
ness H can then be deduced: H = 14.6. Inset of sub-figure (a) 
(resp. inset of sub-figure (b)): plot of the profile of the veloc- 
ity v x (z) (resp. velocity gradient d z v x (z)) in semilogarithmic 
scales. In both insets, the plain straight line corresponds to 
an exponential fit of caracteristic decay length A ~ 3.4. 

as: 

Q = is(z)v x (z)dz (5) 

J z = -R a 

Its variation as well as the one of the flowing layer thick- 
ness H and the mean slope 9 with respect to fl are re- 
ported in Tab|J 

Velocity profiles obtained in the center of the drum 
for various Q are represented in Fig. [SJi. Apart from the 
flowing layer thickness H , the streamwise velocity profile 
at the center of the drum is characterized by two param- 
eters, namely the characteristic decay length A of the ex- 
ponential creep within the static phase and the constant 
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FIG. 4: Velocity profiles obtained for D. = 6 rpm at five dif- 
ferent locations x. The velocity gradient within the flowing 
layer as well as the exponential decreasing within the static 
phase are found to depend weakly on the precise location in 
the drum. They are consequently independent of both the 
flowing layer thickness H and its x derivative d x H 
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TABLE I: Variation of the nondimensionalized flow rate Q, 
the mean angle 9 of the flow and the flowing layer thickness 
H with respect to the rotating speed fi within the elementary 
slice E at the center of the drum. 



shear rate 70 within the flowing layer. Their evolution 
with respect to the flow rate Q is reported in Fig 0J>, c. 
Within the errorbars, A is independent of Q, of the order 
of 3 ± 0.3. This behaviour is similar to what is reported 
in both 3D heap flows exp eriments [52j and 3D rotating 
drum experiments jU H3, , where A was found to be 
A ~ 1.4 and A ~ 2.5 respectively. 

In our numerical simulation, 70 exhibits a weak depen- 
dence with Q (see Fig- Et) - It ranges typically from 0.1 to 
0.25 when Q is made vary from 8 to 58. This dependency 
is compatible with the ones observed experimentally in 
2D rotating drum by Rajchenbach [2!|, who proposed 
that 70 scales as 70 oc (sin# — sin <I>) 1 / 2 / cos 1 / 2 $, where 
4> refers to the Coulomb friction angle. Value of this an- 
gle can be estimated from the variation of the mean flow 
angle with respect to Q (see Fig. Eland next section) and 
was found to be $ = 17.4°. Inset of Fig. shows that 
the scaling proposed by Rajchenbach is compatible with 
our results. It is worth to mention that 70 was found to 
be constant, around 0.5, inde pen dent of Q in 3D experi- 
ments in Hele-Shaw drums ylUJj- This strongly suggests 
some non-trivial effect of either the lateral confinement 
or the flow dimension on the profile within the flowing 
layer. This will be explored in future 3D simulations of 
rotating drums. 
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FIG. 5: (a) Velocity profiles at the center of the drum for 
various rotating speed Q (b) Characteristic decay length A of 
the shear strain d z v x as a function of the flow rate Q. Within 
the errorbars, A is constant, around 3. (c) Constant shear rate 
70 within the flowing layer as a function of the flow rate Q. 
Inset: 7 vs. (sin 9 - sin$) 1/2 / cos ' $ where the Coulomb 
friction angle $ = 17.4° has been identified with the value 
fi e ff(Q = 0) defined in Fig. [(J). The straight line is a linear 
fit 7= 0.75(sin6» - sin $) 1/2 / cos 1/2 $. 



C. Flowing layer thickness and mean flow angle 

The thickness of the flowing layer H is plotted as a 
function of the flow rate Q in Fig. As observed exper- 
imentally urn, 

H scales as \/Q, which is expected since 
the shear rate varies weakly within the flowing layer. 

The mean flow angle 6 can then be assimilated to an ef- 
fective friction coefficient /u e y f = tan 9 between the flow- 
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FIG. 6: (a) Flowing layer thickness H at the center of the 
drum as a function of the flow rate Q. Inset, H vs. s/Q. The 
straight line is a linear fit H = 3\/Q. (b) Variation of the 
effective friction coefficient fi err — tan 9 of the surface flow 
with respect to Q (non-dimensionalized units). The errorbars 
show the standard deviation over the sequence at constant Q. 
The straight line is a linear fit: fi e ff(Q) = 0.31 + 1.9.10" 3 Q. 
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FIG. 7: Mean angular velocity profile ui(z) (averaged over 400 
frames) at the center of the drum obtained for Q — 6 rpm. The 
errorbars correspond to a 99% confident interval. The vertical 
mixed line (resp. the vertical dotted line) show the free surface 
boundary as defined in Fig. (resp. the static/flowing inter- 
face as defined in Fig. Inset, lo vs vorticity V x v = d z v x 
for Q — 2 rpm, (o), fi = 4 rpm (*), £1 = 5 rpm (□), Q — 6 rpm 
(A), Q, — 10 rpm (o), fi = 15 rpm (>). The straight line is 
given by lo — x v. 



particles is equal to half the vorticity. Such relationship 
was observed in Molecular Dynamics simulations of di- 
lute granular flo ws 154. 1 53 , but expected to fail at higher 
volume fraction [&|||&]|. In this latter case, grains were 
expected to organize in layers the grains of which rotate 
in the same direction. This would decrease the mean an- 
gular velocity of the grains, and lu would be smaller than 
x v. In our numerical simulation, such behaviour is 
not observed, which suggests that the grains spins do not 
organize in layer despite the high density of the flow. 



ing layer and the static phase |2J]. Its evolution with 
respect to the flow rate Q is represented in Fig. [fJl The 
effective friction coefficient is found to increase with Q. 
Similar increasing was observed experimentally, - at a 
much larger scale -. It was attributed to wall effects 
since this dependency was observed to be weaker when 
the drum thickness is increased 0, • No such wall ef- 
fects can be invoked in the present study. In other words, 
part of the increase of the effective friction with flow rate 
cannot be induced by wall friction contrary to what was 
suggested in 0, [2j| and should be found in the granular 
flow rheology. 



D. Angular velocity profiles 

A typical mean angular velocity profile u>(z) has been 
represented in Fig. It is interesting to plot lu with re- 
spect to the vorticity V x v = d z v x (see Inset of Fig- EJ - 
There is a clear relationship between these two quanti- 
ties: to = \ V x v in the whole packing independently of Q. 
This relationship is analogue to the one obtained in classi- 
cal hydrodynamics where the mean rotating speed of the 



IV. STRESS ANALYSIS 
A. Stress tensor profile 

Let us now look at stress profiles - that cannot be mea- 
sured experimentally. The stress tensor a is the sum of 
three contributions: a = a c + a k + a r where a c , a k and a r 
refer to the contact, kinetic and rotational components 
of the stress tensor respectively. In our dense free surface 
flows, o~ k and o~ r are found to be negligible with regard 
to a c . One can thus assume that a ~ a c . 

Components of the contact stress tensor associated to 
each bead has been computed for each snapshot of each 
numerical experiment (see Sec[nJ|. The profile of the con- 
tinuum value of each component of the contact stress ten- 
sor - and consequently the total stress tensor - a a p(z) is 
then defined over the elementary slice E located in the 
center of the drum (see Fig. according to the same 
procedure used to calculate velocity and volume frac- 
tion profile. The tensor a is found to be symmetric, i.e. 
o~xz = o~zx- For 2D surface flows, it is thus defined by 
three independent components o~ xx , o~ xz and azz. Typi- 
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cal profile of these components with respect to the depth 
z at the center of the drum are represented in Fig. |2| 

Shapes of these profiles are quite surprising. In the 
rotating frame 5ft, velocity and volume fraction profiles 
were found to be invariant along infinitesimal translation 
along e x within the elementary slice S. In other words, for 
x = 0, one gets dvjdx ~ dv/dx ~ 0. More generally, it is 
commonly assumed that at the center of the drum, the x 
derivative of the stress tensor vanishes 0, Ufl Ll3 ■ 

The Cauchy equations would then read: 



(a) ^ = -j,sin<2 

(b) ^f** = v cos 6 + i>Qv x 



(6) 



where denotes the mean flow angle. The second right- 
handed term of Eq. [BJd is the Coriolis term. This term is 
maximum at the free surface where it reaches 15% of the 
first right-handed gravity term e.g. for = 6 rpm. the 
last right-handed term of Eq. 03 is the centrifugal term. 
This term is maximum at the drum boundary where it 
reaches 1% of the first right-handed gravity term e.g. for 
= 6 rpm. Finally, inertial effects can be neglected and 
the Cauchy equations for pure steady homogenous flows 
would come down to: 



(a) 2g«* = -vsinO 

(b) ^ = -vcos6 



(7) 



and, since the volume fraction v is almost constant, close 
to the random close packing value v RCP = 0.82: 



(a) a xz (z) 

(b) a zz (z) 



_ ZV RCP g j n ( 
_ zl/ RCP cos ( 



(8) 



These predictions were compared to the measured pro- 
files (Fig. |SJ. The measured profile a zz fits well with 
Eq. However, a xz departs from Eq. [SJi within the 
static phase. To understand this discrepancy, one looks 
at the gradient of the stress tensor (see Fig. . The x- 
derivative of the various components were calculated by 
translating the elementary slice S of an increment Sx — 5 
from one side to the other of the reference position x = 0. 
We checked that the obtained values do not depend on 
Sx. Both da zz /dx and da xz /dx vanish within £ at the 
center of the drum. However do~ xx /dx does not. In other 
words, steady surface flows in rotating drums cannot be 
considered as quasi-homogenous even at the center of the 
drum. The Cauchy equations should then read: 



(a) 



da* 



-v cos 9 



(9) 



This may explain the slight discrepancies observed be- 
tween homogenous steady heap surface flows and steady 
surface flows in rotating drum (see e.g. J] for related 
discussions). 



V. CONSTITUTIVE LAWS 
A. Inertial number I 

It was recently suggested 0, H3> EH Eil that the 
shear state of a dense granular flow can be character- 
ized through a dimensionless number /, referred to as 
the inertial number, defined as: 



I = 



d z v x 



(10) 



This parameter can be regarded as the ratio between the 
typical time of deformation l/d z v x and the typical time 
of confinement 1 / y/a zz 

A typical profile of the inertial number I is plotted 
in Fig. HOh ,. This non-dimmensionalized parameter was 
shown to be the relevant parameter to account for the 
transition from the quasi-static regime to the dense in- 
ertial regime in plane shear configurati on, annular shear 
and inclined plane configuration [3, |53, |6(j • Therefore, it 
is natural to consider / as the relevant parameter to de- 
scribe the transition from the quasi-static phase and the 
flowing layer in the surface flow geometry. To check this 
assumption, we determine the value 1th of the inertial 
number at the interface between the static phase/flowing 
layer interface - defined by extrapolating the linear ve- 
locity profile of the flowing phase to zero (see Fig. |3J)- for 
all our numerical experiments carried out at various O. 
Variations of Ith as a function of Q is represented in Fig. 
110b . This threshold is found to be constant, equal to: 



Ii 



Hi 



1.8.10" 



(11) 



which provides a rather strong argument to consider this 
non-dimensionalized parameter as the relevant one to de- 
scribe surface flows. 



B. Rheology 

Now that a relevant parameter describing the local 
shear state of the flow has been proposed, one can dis- 
cuss in more detail the flow rheology. As a first guess, 
it is tempting to consider local constitutive laws relating 
the components of the stress tensor to / through a one- 
to-one relation. In this case, dimensional analysis leads 
to: 



o xz jo zz = fi(I), <7 xx /a zz = k(I) 



(12) 



Typical variations of the effective friction coefficient /j, 
as a function of the inertial number / are plotted on Fig. 
II lb . A semilogarithmic representation (see inset of Fig. 
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FIG. 8: Profiles of the three independent components of the static stress tensor a zz (a), a xx (b) and a xz (c) for £1 = 6 rpm. 
The errorbar correspond to a 99% confident interval. In subfigure (a), (resp. subfigure (c)) the straight line is given by a zz = 
-v RCP cosO (resp. by a xz = — v RCP sin 9) as expected from Cauchy equations for an homogenous surface incompressible flow 
of volume fraction v RCP = 0.82. For Q = 6 rpm, 9 was measured to be 9 = 19.7°. The vertical mixed lines and the vertical 
dotted lines show the free surface boundary and the static/flowing interface as defined in Fig. |5] and Fig. |3] respectively. 




z z z 



FIG. 9: Top: Profiles of the 6 independent component of the contact stress tensor gradient, namely du zz /dz (a), da xx /dz (b), 
da xz /dz (c), da zz /dx (d), da xx /dx (e) and da xz /dx (f) for = 6 rpm. The errorbars correspond to a 99% confident interval. 
In subfigure (a), (c), (d) and (f), the horizontal straight line is given by da zz /dz = — v RCP cos9, da xz /dz — — u RCP sin9, 
da zz /dx = and da xz Jdx = respectively, as expected from Cauchy equations for an homogenous surface incompressible flow 
of volume fraction v = 0.82 and mean flow angle 9 — 19.7° as measured for f2 = 6 rpm. The vertical mixed lines and the 
vertical dotted lines show the free surface boundary and the static/flowing interface as defined in Fig.|5|and Fig. Irrespectively. 



II lb ) shows that data collected for different flow rates Q 
collapse relatively well within the scaling: 



fi = a + b\ogI (13) 

with a ~ 0.35 and b ~ 0.013 when / ranges from 10~ 4 
to 10 _1 . A departure from this scaling is observed when 
I becomes smaller than 10~ 4 . In this latter case, fi de- 
creases more rapidly with /. It is worth to note that the 
scaling given by Ea ll3l is quantitatively similar to the 
one observed in the incline plane geometry [58|, which 
suggests that both free surface flow and flows down to a 
rough incline plane may be described through the same 



constitutive laws. Relating fj, and I through a local con- 
stitutive law seems thus to be relevant. 

Figure ITTb shows the variations of k = cr xx /o~ zz as a 
function of /. In the flowing layer i.e. when / exceed I t h, 
k — > 1. The non monotonic behaviour observed in the 
static phase is much more suprising: The parameter k 
starts from a value lower than 1 at the drum boundary 
k(I — > 0) ~ 0.8, increases and reaches a maximum for I ~ 
10 -3 where k(I ~ 10 -3 ) ~ 1.2 and finally decreases for 
increasing I and tends to 1 within the flowing layer. Such 
observation is very different from the k = 1 behaviour 
observed in the whole materials in both annular shear 
and incline geometry 0, |3^, |58| . 

While the profile {(i(z)} is observed to be invariant 
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FIG. 10: (a): Profiles ol the inertial parameter I(z) for £7 = 
6 rpm. The errorbar correspond to a 99% confident interval. 
The vertical mixed lines and the vertical dotted lines show 
the free surface boundary and the static/flowing interface as 
defined in Fig. |2] and Fig. |3] respectively, (b): Variation of the 
value I t h of the inertial number at the static/flowing interface 
with respect to the flow rate Q. 



along infinitesimal translation, the profile {k(z)} is not 
fFig.ll2|). The ^-derivative of k is found to be almost con- 
stant dk/dx ~ —0.05 within the whole packing. In other 
words, the flow cannot be considered as homogeneous at 
the center of the drum as regard with the parameter k. 
Furthermore, while the curves fj,(I) collected for different 
flow rates Q collapse fairly well, the curves k(I) do not. 
This strongly suggest that the non-local effects implied 
e.g. by the existence of multi-scale rigid clusters embed- 
ded in the flow [33| should be found in the constitutive 
law k(I) rather than in /x(J) contrary to what was sug- 
gested in n h m m . 



FIG. 11: (a): Variation of the effective friction coefficient n = 
<y xz /a zz as a function of the inertial number I for £7 — 6 rpm. 
Inset: Variation of /x as a function of I for £1 = 2 rpm, (o), £7 = 
4 rpm (*), £7 = 5 rpm (□), £7 = 6 rpm (A) and £7 = 10 rpm 
(o) in semilogarithmic scale. The dash-dot line is given by 
/i = 0.35 + 0.013 log I. (b): Variation of the ratio k = a xx /a zz 
as a function of the inertial number I for Q — 6 rpm. Inset: 
Variation of k as a function of / for for Q — 2 rpm, (o), 
£2 = 4 rpm (*), ft = 5 rpm (□), £7 = 6 rpm (A) and £7 = 
10 rpm (o) in semilogarithmic scale. The dash-dot horizontal 
line corresponds to k = 1 
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VI. FLUCTUATION ANALYSIS 

Let us now analyse the fluctuations Sv and 8ui of the 
velocity and the vorticity respectively. Calling c(x, t) the 
"instantaneous" velocity of a bead located at the posi- 
tion x within the elementary slice E at a given time t, 
the fluctuating part of the velocity <5c(x, t) is defined as 
<5c(x, t) = c(x, t) — v x {z)e x where v x {z) denotes the av- 
erage velocity at the depth z (see Fig. [3]}). Profiles of 



-60 -40 -20 o 

z 

FIG. 12: Profiles of the a;-gradient dk/dx of the ratio k = 
Cxx/czz at the center of the drum as observed for £7 = 6 rpm. 
The vertical mixed lines and the vertical dotted lines show 
the free surface boundary and the static/flowing interface as 
defined in Fig. and Fig. 03 respectively. 
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velocity fluctuation Sv 2 (z) are them computed by divid- 
ing £ into layers of one mean bead diameter wide, and 
averaging Sc 2 over all the beads of the 400 frames whose 
center of mass is inside the corresponding layer. Same 
procedure is applied to determine the profiles of angular 
velocities fluctuations. In our athermal granular systems, 
the only time scale is provided by the velocity gradient 
d z v x . Therefore, we looked at profiles of {Sv/d z v x } and 
{Stu/d z v x } rather than direct profiles of {Sv} and {Sco}. 

Figure El displays both translational velocity fluctua- 
tion profile (Fig. 118b) and angular fluctuation (Fig. 113b) 
nondimensionalizcd by the shear rate d z v x . In both cases, 
the nondimensionalizcd fluctuations are found to be con- 
stant within the flowing layer i.e. : 



0; <V 
OUJ 

i>~, t'.r 



2.65 
3.35 



for 
for 



z > —H or 
z > -H or 



I > hh 



(14) 



In the static phase, both 5v/d z v x and 8uj/d z v x are found 
to increase with the distance from the static /flowing in- 
terface. Figure [21 plots both Sv/d z v x (Fig. 114b ) and 
8uj/d z v x (Fig. 114b ) as a function of the inertial num- 
ber /. It evidences the existence of two different scalings 
within the static phase, namely: 
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(15) 



Such scaling are very similar to the one observed in the 
shear geometry [5^ . The importance of these fluctuations 
with regards to the typical rate of deformation d z v x (up 
to 40), as well as the scaling given by Eg. 115b exhibited 
within the static phase are compatible with the picture 
presented in p| to describe quasi-static flow: The aver- 
age grains motion is made of a succession of very slow 
motions when the particle climbs over the next one, and 
a rapid motion when it is pushed back into the next hole 
by the confining picture. 



VII. CONCLUDING DISCUSSION 

Rheologies of 2D dense granular flows were investi- 
gated through Non Smooth Contact Dynamics simula- 
tions of steady surface flows in a rotating drum. Pro- 
files of the different continuum quantities were mea- 
sured at the center of the drum where the flow is 
non-accelerating. Volume fraction v is found to be al- 
most constant, around the Random Close Packing value 
v rcp ^ Q_g2 within the whole packing, except for a tiny 
dilation (few percents) within the flowing layer, as ex- 
pected from dilatancy effects. As observed experimen- 
tally [23, 123, l2a [23, lLLI lM , the streamwise veloc- 
ity profile {v x (z)} is found to be linear within the flow- 
ing layer, and to decrease exponentially within the static 
phase. Mean profile of the angular velocity was also mea- 
sured at the center of the drum and was shown to be 
equal to half of the vorticity in the whole packing. 



FIG. 13: Profile of translational velocity fluctuations Sv (a) 
and angular velocity fluctuation Suj (b) non-dimensionnalized 
by the shear rate d z v x at the center of the drum obtained for 
Q = 6 rpm. The vertical mixed lines and the vertical dotted 
lines show the free surface boundary and the static/flowing 
interface as defined in Fig. [5] and Fig. [3] respectively. 



In a second step, profiles of the three independent com- 
ponent of the stress tensor were measured at the center 
of the drum. Quite surprisingly, the flow is found to be 
non-homogeneous at the center of the drum with regard 
to one of this component, namely o~ xx . In other words, 
d x cr xx does not vanish whereas d x u, c^v, d x o~ zz and d x o~ xz 
vanish. 

The inertial number / - defined as the ratio between 
inertial solicitations and confinement solicitations was de- 
termined. This number is shown to be the relevant one 
to investigate quantitatively the rheology of the surface 
flows. The transition from the static to the flowing phase 
is found to occur to a fixed value Ith of I, independently 
of the flow rate Q. Constitutive laws relating the com- 
ponents of the stress tensor to I were determined. The 
effective friction/i = cr xz /o~ zz is found to increase loga- 
rithmically with /, independently of the flow rate Q. This 
relation is found to match quantitatively the one observed 
in rough incline geometry. On the other hand, the ratio 
k = cr xx /a zz is found to be be significantly different from 
k = 1 in contrast to what was observed in plane shear, 
annular shear, and rough incline geometry [3^, |^. To 
be more precise, k if found to vary non monotonically 
with /. Moreover, d x k is found not to vanish contrary to 
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FIG. 14: (a) Non-dimensionalized velocity fluctuation Sv/d z v x 
as a function of the inertial number I for £1 = 2 rpm, (o), 
£2 = 4 rpm (*), £2 = 5 rpm (□), £7 = 6 rpm (A) and 
£2 = 10 rpm (o). The axes are logarithmic. The slope of the 
plain straight line is —1/2 (b) Non-dimensionalized velocity 
fluctuation Su>/d z function of the inertial number I for 

£7 = 2 rpm, (o), £2 = 4 rpm (*), £2 = 5 rpm (□), £2 = 6 rpm 
(A) and £2 = 10 rpm (o). The axes are logarithmic. The slope 
of the plain straight line is —1/3. In both graphes, the vertical 
dash-dot line show the static/flowing interface as defined by 
/ = hh. 



lation between fi and / would have naturally implied a 
Bagnold velocity profile, as observed in rough incline ge- 
ometry, but not in the present free surface flow geometry. 
In other words, the selection of the velocity profile resides 
more in the function k(I, ...) than in 

Dependencies of {k(I)} with Q, as well as the fact that 
d x k does not vanish lead us to conjecture that the ratio k 
encodes the structure of the percolated network of grains 
in extended contact with each others - referred to as the 
arches network. In this scenario, the structure of this net- 
work - and therefore the ratio k - is related to the global 
geometry of the packing as well as to the orientation of 
the flow. This picture is broadly consistent with nonlocal 
models based on the coexistence of particle chains and 
fluidlike materials However, a more detailed study 

is needed to verify this and understand how k can be 
related to the global structure of the force network. 

Finally, both velocity 5v and angular velocity 5lo 
fluctuations were analysed. These quantities non- 
dimensionalized by the shear rate d z v x were found to 
be constant - independent of the flow rate - within the 
flowing layer thickness. In the static phase, both Sv/d z v x 
and 5v/d z v x were found to decrease as different power- 
laws with I. This behaviour is consistent with the idea 
of an intermittent dynamics generated from the succes- 
sion of rapid rearrangements and slow displacement [j| . 
This change of behaviour at the static / flowing interface is 
broadly consistent with recent measurements of Orpe and 
Khakhar |6l| revealing a sharp transition in the rms ve- 
locity distribution at this interface. Understanding what 
set precisely the scaling laws require precise statistical 
analysis of beads velocities at the grain scale. This rep- 
resents interesting topic for a future investigation. 
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